Descargate este archivo R-Markdown. Este archivo contiene todo el código de esta página.
Para este tutorial necesitaremos los siguientes paquetes. Importante: El paquete RevGadgets esta disponible directamente con install.packages(“RevGadgets”) pero vamos a revisar una versión especial que es la que grafica mapas estocásticos. Por esa razón vamos a bajar la versión que se indica aquí, así que. sigue las instrucciones como se indican abajo.
Revisemos la convergencia básica
# Con el paquete coda
mcmc_run1 <- readTrace(path ="hisse_polinizacion_run_1.log", burnin = 0.1)
## Reading in log file 1
mcmc_trace <- as.mcmc(mcmc_run1[[1]])
traceplot(mcmc_trace)
effectiveSize(mcmc_trace)
## Iteration Posterior Likelihood
## 0.0000 1570.0326 4876.6295
## Prior extinction[1] extinction[2]
## 1551.4522 555.6097 439.0926
## extinction[3] extinction[4] extinction_alpha[1]
## 665.4208 712.8727 484.9212
## extinction_alpha[2] extinction_beta[1] extinction_beta[2]
## 305.5809 360.2756 234.7382
## net_diversification[1] net_diversification[2] net_diversification[3]
## 2265.9664 637.0799 6278.9435
## net_diversification[4] q_0A0B q_0A1A
## 5897.2707 7155.9250 4266.2468
## q_0B0A q_0B1B q_1A0A
## 3418.0186 18113.4734 3455.6445
## q_1A1B q_1B0B q_1B1A
## 3693.0056 5912.6036 3784.2706
## root_frequencies[1] root_frequencies[2] root_frequencies[3]
## 24185.6143 73529.7285 72125.3665
## root_frequencies[4] speciation[1] speciation[2]
## 71765.3892 804.1174 1036.4692
## speciation[3] speciation[4] speciation_alpha[1]
## 3835.6708 4660.1254 853.2841
## speciation_alpha[2] speciation_beta[1] speciation_beta[2]
## 1076.8263 877.0545 472.4078
Ahora hagamos un mejor trabajo con los gráficos de convergencia
# Esta parte se hace con ggplot2 y con RevGadgets
mcmc_run1 <- readTrace(path ="hisse_polinizacion_run_1.log", burnin = 0.1)
## Reading in log file 1
mcmc_run1<-data.frame(mcmc_run1)
mcmc_run1<- cbind(mcmc_run1,run=rep("run 1",length(mcmc_run1$Iteration)))
mcmc_run2 <- readTrace(path ="hisse_polinizacion_run_2.log", burnin = 0.1)
## Reading in log file 1
mcmc_run2<-data.frame(mcmc_run2)
mcmc_run2<- cbind(mcmc_run2,run=rep("run 2",length(mcmc_run2$Iteration)))
mcmc_table<-rbind(mcmc_run1,mcmc_run2)
trace_plot<- ggplot(mcmc_table, aes(x=Iteration,y=Posterior,group=run))+
geom_line(aes(color=run))+
theme_classic()
trace_plot
#Aqui se nota que hay que cortar un poco mas de burn-in (15,000 generaciones), asegúrense siempre de cortar mas valores si hace falta
Queremos graficar y obtener estadísticas resumen de los parametros de interés. Recordemos que para el modelo hisse los parámetros que estimamos van a ser los siguientes:
Tasas de transición\[q_{0A1A}, q_{1A0A}, q_{0B1B}, q_{1B0B},q_{0A0B}, q_{0B0A}, q_{1A1B}, q_{1B1A}\]
Tasas de transición entre estados escondidos \(q_AB, q_BA\) (se llaman hidden_rate en el código de RevBayes)
Tasas de diversificación neta \[r_{0A}=\lambda_{0A}-\mu_{0A},r_{1A}=\lambda_{1A}-\mu_{1A},r_{0B}=\lambda_{0B}-\mu_{0B}, r_{1B}=\lambda_{1B}-\mu_{1B}\]
Alternativo: la fracción de extinción \[\epsilon_{0A}=\mu_{0A}/\lambda_{0A}, \epsilon_{1A}=\mu_{1A}/\lambda_{1A}, \epsilon_{0B}=\mu_{0B}/\lambda_{0B},\epsilon_{1B}=\mu_{1B}/\lambda_{1B}\]
Frecuencias de la raíz
Para los modelos más complicados siempre prefiero graficar violines, y tener más control en los colores y la figura, por eso mi selección de ggplot2. La comparación de violines es super útil.
# Esta parte se hace con ggplot2
traitcols <- c("#3D348B", "#5A4FCF","#F18701", "#F35B04", "#1B7F3A", "#4CAF72","#D4A300", "#F2D15C")
hisse<- read.table("hisse_polinizacion_run_1.log", header=TRUE)
hisse<- hisse[-seq(1,15000,1),] # nunca olvidemos quitar el burn-in con esta instrucción tengo mas control de cuantas iteraciones remuevo.
tasas_transicion <- data.frame(dens=c(
hisse$q_0A1A,hisse$q_1A0A,
hisse$q_0B1B,hisse$q_1B0B,
hisse$q_0A0B,hisse$q_0B0A,
hisse$q_1A1B,hisse$q_1B1A),
rate=rep(c("q_0A1A","q_1A0A","q_0B1B","q_1B0B","q_0A0B","q_0B0A","q_1A1B","q_1B1A"),each=length(hisse$q_0A1A))
)
tasas_transicion$rate <-
factor(tasas_transicion$rate,
levels = c("q_0A1A","q_1A0A","q_0B1B","q_1B0B","q_0A0B","q_0B0A","q_1A1B","q_1B1A"))
violin_transitions<- ggplot(tasas_transicion,aes(x=rate,y=dens, fill=rate))+
geom_violin(trim=FALSE)+
labs(title="Tasas de transicion")+
scale_fill_manual( values = traitcols)+
theme_classic()
violin_transitions
¿Cuál es la interpretación del gráfico?¿Existen diferencias? R:Tu Respuesta
Ahora vamos comparar Transiciones entre estados escondidos
## Vamos a crear una tabla de las diferencias entre las transiciones
test_tasas_transicion_hisse<-
data.frame(
dens=c((hisse$q_0A1A - hisse$q_1A0A),(hisse$q_0B1B - hisse$q_1B0B),
(hisse$q_0A1A - hisse$q_0B1B),(hisse$q_1A0A - hisse$q_1B0B),
(hisse$q_0A0B - hisse$q_0B0A),(hisse$q_1A1B - hisse$q_1B1A),
(hisse$q_0A0B - hisse$q_1A1B),(hisse$q_0B0A - hisse$q_1B1A)),
difference=rep(c("q_0A1A - q_1A0A","q_0B1B - q_1B0B",
"q_0A1A - q_0B1B","q_1A0A - q_1B0B",
"q_0A0B - q_0B0A","q_1A1B - q_1B1A",
"q_0A0B - q_1A1B","q_0B0A - q_1B1A"
),each=length(hisse$q_0A1A)))
test_tasas_transicion_hisse$difference <-
factor(test_tasas_transicion_hisse$difference,
levels = c("q_0A1A - q_1A0A","q_0B1B - q_1B0B",
"q_0A1A - q_0B1B","q_1A0A - q_1B0B",
"q_0A0B - q_0B0A","q_1A1B - q_1B1A",
"q_0A0B - q_1A1B","q_0B0A - q_1B1A"))
violin_hidden <- ggplot(test_tasas_transicion_hisse,aes(x=difference,y=dens, fill=difference))+
geom_violin(trim=FALSE,draw_quantiles = c(0.025, 0.975))+ #Intervalos de maxima probailidad
labs(title="Tasas de transición entre estados escondidos ")+
geom_hline(yintercept = 0, linetype = "dashed", lwd = 1)+
scale_fill_manual( values = traitcols)+
theme_classic()+
theme(legend.position = "none", plot.title = element_text(hjust = 0.5),axis.text.x = element_text(angle = 45, hjust = 1))
## Warning: The `draw_quantiles` argument of `geom_violin()` is deprecated as of ggplot2
## 4.0.0.
## ℹ Please use the `quantiles.linetype` argument instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
violin_hidden
¿Cuál es la interpretación del gráfico?¿Existen diferencias? R:Tu Respuesta
# Esta parte se hace con ggplot2
divcols<-c("#E63946","#1D3557","#F3A5AB", "#457B9D")
# Recuerda de RevBayes 1=0A, 2=1A, 3=0B, and 4=1B
tasas_diversificacion_neta <- data.frame(dens=c(hisse$speciation.1.-hisse$extinction.1.,hisse$speciation.2.-hisse$extinction.2., hisse$speciation.3.-hisse$extinction.3.,hisse$speciation.4.-hisse$extinction.4.) ,rate=rep(c("net_div_0A","net_div_1A","net_div_0B","net_div_1B"),each=length(hisse$speciation.1.)))
tasas_diversificacion_neta$rate<- factor(tasas_diversificacion_neta$rate,levels=c("net_div_0A","net_div_1A","net_div_0B","net_div_1B"))
violin_diversificacion<- ggplot(tasas_diversificacion_neta,aes(x=rate,y=dens, fill=rate))+
geom_violin(trim=FALSE)+
labs(title="Tasas de diversificación neta")+
scale_fill_manual( values = divcols)+
theme_classic()
violin_diversificacion
¿Qué podemos concluir de este gráfico?
Una de las ventajas más importante de las distribuciones posteriores es que rápidamente podemos concluir si dos parámetros son iguales o distintos. Formalicemos entonces la prueba de hipotésis para saber si el caracter esta asociado a la diversificación.
Construimos dos estadísticas de resumen que son las diferencias entre la diversificación neta de 0 y 1 para el estado A y el estado B respectivamente:
Si estas diferencias son 0 con alta probabilidad esto significa que la polinización y sus estados no hacen la diferencia en el proceso de diversificación.
# Esta parte se hace con ggplot2
difcols<-c("#FF006E","#FFC2DC")
T_diferencias<- data.frame(dens=c((hisse$speciation.1.-hisse$extinction.1.)-(hisse$speciation.2.-hisse$extinction.2.),(hisse$speciation.3.-hisse$extinction.3.)-(hisse$speciation.4.-hisse$extinction.4.)),diferencia=rep(c("T_A","T_B"),each=length(hisse$speciation.1.)))
violin_difference<- ggplot(T_diferencias,aes(x=diferencia,y=dens, fill=diferencia))+
geom_violin(trim=FALSE,draw_quantiles = c(0.025, 0.975))+
labs(title="Estadisticas resumen")+
scale_fill_manual( values = difcols)+
geom_hline(yintercept = 0,linetype="dashed",lwd=1)+
theme_classic()
# Observemos que 0 cruza la diferencia entre las tasas lo que quiere decir que pueden ser iguales
violin_difference
¿Cuál es la interpretación del gráfico?¿Existen diferencias? R:Tu Respuesta
Formalizando la prueba de hipótesis \[ H_0: T_A=0 \textrm{ y } T_B=0\]
# Esta parte se hace con dplyr
cuantiles <- T_diferencias %>% group_by(diferencia)%>%reframe(res=quantile(dens,probs=c(0.025,0.975)))
cuantiles
## # A tibble: 4 × 2
## diferencia res
## <chr> <dbl>
## 1 T_A -0.284
## 2 T_A 0.390
## 3 T_B -0.636
## 4 T_B -0.0135
Observamos que el intervalo de credibilidad del 95% de la distribución de \[T_A\] es (-0.284, 0.389) y el intervalo de credibilidad del 95% de la distribución de \[T_B\] es (-0.636, -0.013). Cómo 0 pertenece al intervalo \[T_A\], significa que el tempo de la diversificación para A es igual, mientras 0 no pertenece al intervalo \[T_B\]entonces significa que el tempo de la diversificación para el estado 1 es diferente. Importante: Noten que en las pruebas de hipótesis bayesiana no utilizo, p-valores, ni significancia, ni rechazo, solo probabilidad e intervalos de credibilidad. Esto sucede porque la manera en que interpretamos probabilidad en bayesiana es distinto. Tengan cuidado cuando interpreten.
Segunda parte: ¿Qué pasa con los estados escondidos?
Construimos dos estadísticas de resumen que son las diferencias entre la diversificación neta de A y B para el estado 0 y el estado 1 respectivamente:
Si estas diferencias son 0 con alta probabilidad esto significa que hay cambios en el tempo de la diversificación pero se deben a algo más que nosotros no medimos en la filogenia.
# Esta parte se hace con ggplot2
difcols<-c("#1DB32C","#BFEEC3")
T_diferencias<- data.frame(dens=c((hisse$speciation.1.-hisse$extinction.1.)-(hisse$speciation.3.-hisse$extinction.3.),(hisse$speciation.2.-hisse$extinction.2.)-(hisse$speciation.4.-hisse$extinction.4.)),diferencia=rep(c("T_0","T_1"),each=length(hisse$speciation.1.)))
violin_difference<- ggplot(T_diferencias,aes(x=diferencia,y=dens, fill=diferencia))+
geom_violin(trim=FALSE,draw_quantiles = c(0.025, 0.975))+
labs(title="Estadisticas resumen")+
scale_fill_manual( values = difcols)+
geom_hline(yintercept = 0,linetype="dashed",lwd=1)+
theme_classic()
violin_difference
¿Cuál es la interpretación del gráfico?¿Existen diferencias? R:Tu Respuesta
Formalizando la prueba de hipótesis \[ H_0: T_0=0 \textrm{ y } T_1=0\]
# Esta parte se hace con dplyr
cuantiles <- T_diferencias %>% group_by(diferencia)%>%reframe(res=quantile(dens,probs=c(0.025,0.975)))
cuantiles
## # A tibble: 4 × 2
## diferencia res
## <chr> <dbl>
## 1 T_0 -0.453
## 2 T_0 -0.195
## 3 T_1 -1.02
## 4 T_1 -0.339
Observamos que el intervalo de credibilidad del 95% de la distribución de \[T_0\] es (-0.452, -0.195) y el intervalo de credibilidad del 95% de la distribución de \[T_1\] es (-1.024, -0.339). Cómo 0 no pertenece a estos intervalo entonces \[P(H_0\lvert Datos)<0.05\] lo que significa que el tempo de la diversificación para A y B es diferente. Este resultado unido al resultado anterior indican que el tempo cambia debido a algo más que nosotros no medimos pero se manifiesta en el árbol. Estos resultados son equivalentes ajustar el modelo de caracter independiente CID y seleccionarlo, indicando que la diversificación es independiente del caracter de interés.
Por que la curiosidad mató al gato chequemos que paso con BiSSE
Baja los siguientes archivos:
# Esta parte se hace con ggplot2
bisse<- read.table("bisse_polinizacion_run_1.log", header=TRUE)
bisse<- bisse[-seq(1,15000,1),]
divcols<-c("#E63946","#1D3557")
# Recuerda de RevBayes 1=0A, 2=1A, 3=0B, and 4=1B
netdiversification_rates<- data.frame(dens=c(bisse$net_diversification.1., bisse$net_diversification.2.) ,rate=rep(c("net_div_0","net_div_1"),each=length(bisse$net_diversification.1.)))
violin_diversification<- ggplot(netdiversification_rates,aes(x=rate,y=dens, fill=rate))+
geom_violin(trim=FALSE)+
labs(title="Tasas de diversificación neta")+
scale_fill_manual( values = divcols)+
theme_classic()
violin_diversification
Prueba estadística
difcols<-c("#1DB32C","#BFEEC3")
T_diferencia<- data.frame(dens=(bisse$net_diversification.1.-bisse$net_diversification.2.),diferencia=rep("T",each=length(bisse$net_diversification.1.)))
violin_difference<- ggplot(T_diferencia,aes(x=diferencia,y=dens, fill=diferencia))+
geom_violin(trim=FALSE,draw_quantiles = c(0.025, 0.975))+
labs(title="Estadìstica resumen")+
scale_fill_manual( values = difcols)+
geom_hline(yintercept = 0,linetype="dashed",lwd=1)+
theme_classic()
# Observemos que 0 cruza en una probabilidad muy pequeña de la distribución la diferencia entre las tasas lo que quiere decir que pueden ser iguales
violin_difference
## 0 no esta en el intervalo de probabilidad del 95%
quantile ((bisse$net_diversification.1.-bisse$net_diversification.2.),probs=c(0.025,0.975))
## 2.5% 97.5%
## -0.558164 -0.092405
¿Qué hubieras concluido para BiSSE? ¿Cuál es la mayor falla de este modelo?